# Required packages
require(tidyverse)
require(survival)
require(tidymodels)
require(flextable)

# Version used:
# survival_3.3-1

# Load data
load("01_Data/spatial_dimension_data.Rdata")

#### 1 -  Modify censored observations data ####
model_survival_0 <- model_base %>% 
  group_by(id) %>% 
  arrange(desc(date)) %>% 
  slice(1)


# Set censored status
model_base$status <- 1
model_survival_0$status <- 0

# Calculate censored interval
model_survival_0$interval <- as.Date("2020-12-31")-model_survival_0$Date 
model_survival_0$interval <- as.numeric(model_survival_0$interval)

# Filter very long intervals
model_survival_0 <- model_survival_0 %>% 
  filter(interval < 1500)

# Combine dataset for model
model_survival <- rbind(model_base, model_survival_0)

#### 2 - Estimate model ####
m1 <-
  (
    coxph(
      Surv(interval, status) ~ risk_factor + time,
      model_survival,
      robust=T
    )
  )

#### 3 - Table A5: Regression table survival model ####
m1_tidy <- tidy(m1)
m1_glance <- glance(m1) %>%
  select(nobs,
         r.squared,
         AIC,
         BIC) %>%
  pivot_longer(cols = everything(),
               names_to = "term",
               values_to = "estimate")

m1_table <- bind_rows(m1_tidy, m1_glance)


m1_table$term[m1_table$term == "time"] <- "Distance (hours)"
m1_table$term[m1_table$term == "risk_factorhigh"] <- "Risk (high)"
m1_table$term[m1_table$term == "risk_factormedium"] <- "Risk (medium)"
m1_table$term[m1_table$term == "r.squared"] <- "(pseudo-)R2"
m1_table$term[m1_table$term == "nobs"] <- "N"

m1_table$stars <- ""

m1_table$stars[(m1_table$p.value) < 0.1] <- "*"
m1_table$stars[(m1_table$p.value) <  0.05] <- "**"
m1_table$stars[(m1_table$p.value) < 0.01] <- "***"

m1_table$estimate <- ifelse(
  !is.na(m1_table$std.error),
  paste0(
    format(
      round(m1_table$estimate, 2),
      nsamll = 2,
      trim = T
    ),
    m1_table$stars,
    " ",
    "(",
    format(
      round(m1_table$std.error, 2),
      nsamll = 2,
      trim = T
    ),
    ")"
  ),
  format(
    round(m1_table$estimate, 2),
    nsmall = 2,
    trim = T
  )
)

summary(m1)
m1_table$estimate[m1_table$term=="N"] <- 2281

m1_table <- m1_table %>% 
  select(term,estimate) 


m1_table <-
  rbind(
    c(
      "DV", "Interval"), m1_table)

m1_table %>% flextable() %>% 
  hline(i = 1, j = -1, border = NULL, part = "body") %>% 
  hline(i = 4, j = -1, border = NULL, part = "body") %>% 
  
  add_footer_lines(" *** p < 0.01; ** p < 0.05; * p < 0.1.
  Robust standard errors in parentheses
  NOTE: This model includes censored observations") %>% 
  fontsize(size=7, part="all") %>% 
  width(j = -1, 3, unit = "cm") %>% 
  width(j = 1, 3, unit = "cm") %>% 
  align(
    j = -1,
    align = "center",
    part = "all",
    ) %>% 
  save_as_docx(path="03_Tables/appendix_A5.docx")






rm(list = ls())
